Stability of Boolean and continuous dynamics 
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Regulatory dynamics in biology is often described by continuous rate equations for continuously 
varying chemical concentrations. Binary discretization of state space and time leads to Boolean 
dynamics. In the latter, the dynamics has been called unstable if flip perturbations lead to damage 
spreading. Here we find that this stability classification strongly differs from the stability properties 
of the original continuous dynamics under small perturbations of the state vector. In particular, 
random networks of nodes with large sensitivity yield stable dynamics under small perturbations. 

PACS numbers: 89.75.Hc,05.45.-a,87.10.-e,45.05.+x 



The functioning of organisms on the molecular level 
is a research topic of increasing attention. Survival and 
reproduction requires an autonomous regulation of chem- 
ical concentrations in the living cell. Modeling such regu- 
latory dynamics, various mathematical approaches have 
been studied, from discrete to continuous methods, from 
deterministic to stochastic techniques, from static to dy- 
namical models, from detailed to coarse grained perspec- 
tives [H, see ref. Q for an overview. 

Boolean dynamics is a framework for modeling 

regulatory systems, especially for precise sequence con- 
trol as observed in morphogenesis @ and cell cycle dy- 
namics [l(| but also in the regulation of the metabolism 
[Tlj . Using binary (on/off) concentrations as an idealiza- 
tion, Boolean dynamics directly implements the logical 
skeleton of regulation. Values of system parameters such 
as binding constants, production and degradation rates 
etc. are not needed. This abstraction simplifies com- 
putation and analytical treatment. Boolean networks 
have been extracted directly from the literature jf| HH 
of known biochemical interactions or obtained by dis- 
cretization of differential equation models (l3| . Known 
state sequences and responses of several systems have 
been faithfully reproduced by the discrete models 0, ■ 

Despite these benefits, modelers do not employ 
Boolean dynamics as widely as ordinary or delay differ- 
ential equations. The latter arc embedded in an estab- 
lished framework for state- continuous dynamical systems 
flij which itself builds on the mathematical foundations 
of linear algebra and infinitesimal calculus. In particular, 
the definition of stability of a solution under small per- 
turbations is based on the consideration of infinitesimally 
small neighborhoods in state space. Stability checks for 
solutions of the dynamical equations are a salient part 
of mathematical modeling. Unstable solutions are not 
expected to be observed in a real-world system. 

In the state- discrete Boolean dynamics, large pertur- 
bations are normally implemented as a flip, where the 
state of a single Boolean variable is inverted. Then the 
evolution of the damage is tracked. The damage is the 



difference between the state of the perturbed and the un- 
perturbed system. The return map of the expected size 
of the damage is known as Derrida plot [15| . Numer- 
ous studies have elucidated the effect of flip perturba- 
tions on regulatory dynamics with Boolean states fill - 
ip . When asking if a gene-regulatory system reproduces 
a prescribed trajectory despite noise, large perturbations 
are to be considered in the case of low copy numbers 
of regulatory molecules and bursty stochastic response 
[22l ]. Small perturbations, however, are more appropri- 
ate when modeling systems with large copy numbers and 
an integrative response to filter out bursts, see e.g. (23j . 

Here we find that the clear distinction between the two 
types of perturbations is crucial. In a continuous system, 
stability or instability under small perturbations is not 
indicative of the effect of flip perturbations. Likewise, 
probing a Boolean system with flip perturbations does 
not necessarily provide information about the stability of 
the continuous counterpart under small perturbations. 

An n-dimensional Boolean map / : {0, 1}" — > {0, l} ra 
gives rise to a time-discrete dynamics 



x(t + l) = f(x(t)) 



(1) 



with x = (xi, . . . ,x n ) being a Boolean state vector (bit 
string) of n entries. Such a map is equivalent to a Boolean 
network. When / is pictured as a network, a node corre- 
sponds to a coordinate i of the Boolean state vector and 
a directed edge j —¥ i (from node j to node i) is present 
if the Boolean function fi explicitly depends on the j-th 
coordinate. 

Let us now define a continuous dynamics whose dis- 
cretization readily leads to the Boolean map in Eq. ((T|). 
Taking values yi(t) € [0,1], i E {1, . . . , n}, t G R, the 
states evolve according to the delay differential equation 



W(i + 1) = a Bgn(/(y(f)) - y<(t + 1)) 



(2) 
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with a an inverse time constant. For large a, this is 
essentially Boolean dynamics with fast but continuous 
switching between the saturation values. The simplest 
choice is / = / o G with the component- wise step func- 
tion, ®i(y) = 1 if yi > 1/2 and &i{y) = otherwise. 
This choice of continuous dynamics is in close correspon- 
dence with the discrete dynamics in the following sense. 
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Suppose x(0), x(l), x(2), ... is a solution of Eq. (QJ. Let 
y(t) be a solution of Eq. © such that there is a time 
interval [ti,^] with y(s) — x(0) for all s £ [ii,i2]- Then 
for all future times t £ N and all s £ [ti, ta] 



s(t) = j/(/3t + s) 



(3) 



with /3 = 1 + l/(2a). The closest resemblance be- 
tween Boolean and continuous dynamics is obtained 
when choosing the same initial condition, that is y(s) = 
x(0) for all s £ [—1,0]. Similar correspondence between 
Boolean maps and ordinary differential equations has 



been studied earlier neglecting transmission delay 24 1 
or implementing more complicated differential equations 
[25l - l28j compared to Equation @. 

Perturbations. - Given a map /, the evolution of 
states is uniquely determined by Eq. ([2]) by an initial 
condition y(t) on a time interval of unit length, here taken 
as [—1, 0] =: /. We restrict ourselves to initial conditions 
that do not vary on /, y(t) = y(0) for all t £ /. An initial 
condition with a small perturbation is generated as 



y[{t) : =e i (l-y i (t)) + (l-e i )y i (t) 



(4) 



for t £ /. The perturbation amplitudes are arbitrary 
numbers £ ]0, l/2[. An initial condition with a flip 
perturbation is generated as 



y\{t) 



l — yi(t) i£i = l 
yi (t) otherwise 



(5) 



for t £ / and an arbitrary node I £ {1, . . . ,n}. Note that 
the total amplitude J2i e « °f a small perturbation may ex- 
ceed the unit amplitude of a flip perturbation. A small 
perturbation produces small deviations from the origi- 
nal state potentially at each node. A flip perturbation 
concentrates a maximal deviation at a single node. 

We say that the system heals from the perturbation 
if the dynamics from perturbed and unperturbed initial 
condition eventually become the same except for an ar- 
bitrary time lag. Formally, healing from a small pertur- 
bation means that there are to > and t > —to such 
that 



y(t)=y'(t + T) 



(G) 



for all t > t . Healing from a flip perturbation means 
that Eq. © holds analogously for y- instead of y' . We 
define the heal time iheai as the smallest time to for which 
this holds. 

Fixed points and bistable circuits. — Let us first con- 
sider a fixed point as the simplest dynamical behaviour. 
A fixed point of the continuous dynamics is a state vector 
y* such that constant y(t) = y* is a solution of Eq. ©. 
This in turn means that the time derivative vanishes at 
all times, equivalent to y* = f(y*). The fixed points of 
the continuous dynamics are exactly the fixed points of 
the discrete map /. A small perturbation to a fixed point 
y* always heals, because values after applying the tresh- 
old 9 remain unchanged, f(y'(t)) = y* for all t £ /. All 




FIG. 1: (color online). Dynamics of two mutually activating 
nodes, (a) State space of the Boolean system described by 
Eq. ([7]). Thin arrows indicate the mapping / of states by 
the dynamics, thick bidirectional arrows stand for flip pertur- 
bations. Indicated by shaded areas, the system has three 
dynamical modes (attractors): two fixed points (0,0) and 
(1, 1) and a cycle of length 2 involving the states (0, 1) and 
(1,0). (b) Time evolution of the corresponding continuous 
system in Equation ([S| with initial condition £i(0) = 1 (thick 
curve) and X2(0) — (thin curve). The two nodes switch in 
a synchronous mode as indicated by vertical double arrows 
akin to the Boolean state sequence (0, 1), (1, 0), (0, 1), ... . (c) 
Time evolution from perturbed initial condition, a;i(0) < 1, 
£2(0) > 0. The perturbation translates into a phase lag in 
switching that does not heal out. 



fixed points are stable under small perturbations. How- 
ever, a flip perturbation to a fixed point does not always 
heal. The bistable switch is an example. Consider a two- 
dimensional map / with f(x\,X2) = {x2,xi)- It gives 
rise to the dynamics 



SCi(t + l) =x 2 (t) X 2 (t + 1) = Xl(t) 



(7) 



with fixed points (0,0) and (1,1). After perturbing a 
fixed point by flipping one node's state, the system does 
not return to the fixed point. It remains in the set of 
the state vectors (0, 1) and (1,0) constituting a limit cy- 
cle, cf. Figure [TJa). The stability of the fixed points is 
not obtained when probing the dynamics with flip per- 
turbations. The bistable switch constitutes a first simple 
example of systems with different stability properties un- 
der flip and small perturbations. 

In the continuous counterpart of the alternating 
Boolean state (0,1) and (1,0), small perturbations do 
not heal, see Figure [1Tb, c). The effect of a small pertur- 
bation is to induce a phase lag in the oscillation, being 
discussed in earlier work [25|, [29l - [3l| . 

Stability in random networks. — We now compare the 
effects of the two types of perturbations on dynamics in 
randomly generated networks. An ensemble of random 
Boolean networks (RBN) [f| is defined by the number 
of nodes n, the number of inputs K of each node, and 
the probability distribution of Boolean functions 7r(/). 
The latter is taken as a maximum entropy ensemble 
7Ta(/) oc exp(As(/)) under a given average sensitivity (s). 
The sensitivity s(f) of a Boolean function / is the num- 
ber of flips at one of the K inputs that lead to a change 
of the output value, averaged over all input vectors [32|. 
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average sensitivity <s> 



FIG. 2: (color online). Stability of dynamics in random net- 
works under perturbation by spin flip (dashed) and under 
continuous perturbation (solid lines) in random networks with 
K — 2 and K — 4 inputs per node. Symbols distinguish sys- 
tem size n = 300 (o), 1000 (□) and 3000 (o). Each data 
point gives the relative frequency of healed out perturbations 
on a set of 10 4 independent random realizations of network, 
initial condition and perturbation. Each amplitude e, of a 
small perturbation is drawn independently from the uniform 
distribution on an interval [0; r] with < r < 0.5. The results 
are independent of the choice of r. As a general invariance of 
the dynamics of Equation with / = / o G, the qualitative 
effect (healing or spreading) of a small perturbation is not al- 
tered when the amplitude vector is multiplied with a positive 
scalar keeping each amplitude e< < 0.5. 



The resulting value s(f) lies in the range from zero (for 
a constant function /) to K, obtained for a parity func- 
tion where for all input vectors, a flip of a single input 
state flips the output. For RBN, where the K inputs of 
each node are drawn randomly and independently from 
the set of n nodes, the average sensitivity (s) is the cru- 
cial parameter determining the system's response to flip 
perturbations [32j|. In the limit n — > oo, these perturba- 
tions heal in ensembles with (s) < 1; they spread when 
(s) > 1. This change of behaviour in dependence of (s) is 
reproduced in Figure [5] (dashed lines) for varying K and 
n. 

As our main result, we show in Fig. [2] that the (s)- 
dependence of the healing probability under flip pertur- 
bations is qualitatively different from that under small 
perturbations. Only in the so-called critical region of 




system size 



FIG. 3: (color online). The average time to heal from a small 
perturbation increases linearly with the number of nodes in 
the system for sensitivity (s) > 1, and sublinearly otherwise. 
The dashed line has slope 1 in this double-logarithmic plot. 
Each data point is the average over theai for the subset of 
healing realizations. Realizations of network, initial condition 
and perturbation are the same as in Figure [21 



(s) w 1, small perturbations spread. Both for (s) <C 1 
and (s) ^> 1, the healing probability tends towards 1. 
This effect is enhanced by increasing system size. In the 
limit of n — > oo one may expect a finite probability of 
non-healing only at (s) = 1. Then the dynamics is al- 
most always stable under small perturbations. 

The average time theai to heal from small perturba- 
tions increases moderately with system size as shown in 
Figure [3J For average sensitivity above 1, we observe a 
linear increase (thoai) oc n. For lower values of the average 
sensitivity, the increase is sublinear. 

The dynamics we have studied so far is simple but not 
the only possibility to pass from the Boolean map to a 
continuous flow. In order to check to what extent our 
results depend on this choice we repeat simulations for 
K = 2 with an alternative function / (cf. Equation @) 
now taking into account cooperative effects between in- 
puts. Figure [4] shows that the same qualitative result 
obtains under this choice, see figure caption for details. 

In summary, we have shown that the dynamics of large 
random networks of switch-like elements typically recov- 
ers from small perturbations of the state vector. Healing 
is observed naturally at low sensitivity. However, also 
large sensitivities of the nodes' functions render the long- 
term behaviour of the whole system insensitive to small 
perturbations. Instability is observed only in an interme- 
diate sensitivity regime that shrinks as systems become 
larger. 

The behaviour under small perturbations is essentially 
different from the established stability diagram for RBN. 
Under flip perturbations, RBN display a transition from 
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average sensitivity <s> 



FIG. 4: (color online). Healing probabilities remain qual- 
itatively the same (cf. Figure [2} when using the alternative 
transfer function fi(y) = Q{hi(y)) with hi(y) = ayjyk+biyj + 
foj/fc + c; for node i taking inputs from nodes j and k. The 
parameters a,b\,b2,c are chosen such that hi(y) = fi(y) for 
J/jiJ/fc S {0, 1}. If, for instance, fi is an AND then a = 1 and 
6i = 62 = c = so = 1 if and only if the product of 

inputs j/jj/fe > 1/2. Each data point is the healing fraction 
of 1000 realizations of given average sensitivity and system 
size n — 30 (o), 100 (□) and 300 (o). The perturbation am- 
plitude €i is drawn from the uniform distribution on [0; 0.01] 
independently for each node i. 



healing to non-healing (damage spreading) behaviour at 
average sensitivity 1. It has been suggested that net- 
works of regulatory switches position themselves at this 
transition (33[, known as the edge of chaos J34[- Then 
some but not all flip perturbations spread and therefore 
allow for complex information processing without render- 
ing the system unreliable under noise. 

According to our findings, a complementary scenario 
is worth discussing. The apparent conflict between re- 
sponsiveness to external input signals and resilience to 
intrinsic noise dissolves when these influences act as per- 
turbations at separate scales: noise corresponds to small 
perturbations whilst input signals are interpreted as the 
flipping of a state. Under these assumptions, noise re- 
silience and responsiveness are compatible rather than 
conflicting in the regime of average sensitivity above 1. 
Systems that combine both beneficial properties are ob- 
tained "for free" in random networks of sufficiently sen- 
sitive switching elements. 
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